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Abstract 

Background: Altered DNA methylation patterns represent an attractive mechanism for understanding the phenotypic 
changes associated with human aging. Several studies have described global and complex age-related methylation 
changes, but their structural and functional significance has remained largely unclear. 

Results: We have used transcriptome sequencing to characterize age-related gene expression changes in the human 
epidermis. The results revealed a significant set of 75 differentially expressed genes with a strong functional relationship 
to skin homeostasis. We then used whole-genome bisulfite sequencing to identify age-related methylation changes at 
single-base resolution. Data analysis revealed no global aberrations, but rather highly localized methylation changes, 
particularly in promoter and enhancer regions that were associated with altered transcriptional activity. 

Conclusions: Our results suggest that the core developmental program of human skin is stably maintained through 
the aging process and that aging is associated with a limited destabilization of the epigenome at gene regulatory 
elements. 
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Background 

Epigenetic mechanisms regulate the interpretation of 
genetic information and are thus intricately linked to 
cellular differentiation and tissue specification. Epigenetic 
mechanisms include DNA methylation and covalent 
histone modifications [1,2]. These mechanisms work in a 
coordinated manner during development and differenti- 
ation to execute specific gene expression programs [3,4]. 

DNA methylation is a conserved epigenetic mechanism 
with a well-known role in cell fate specification [5,6]. In 
mice, this function appears to be essential for organismal 
development, as genetic deficiencies for DNA methyl- 
transferases cause severe developmental phenotypes 
and embryonic lethality [7,8]. More recently, significant 
attention has also been focused on adaptive functions 
of DNA methylation, which could provide the foundation 
for the integration of environmental signals [9]. Initially 
observed as epigenetic variations between monozygotic 
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twins [10], age-related methylation differences have 
now been described in several independent studies and 
tissues [11]. The results suggested that aging induces 
global and complex changes in the human methylation 
landscape [12-15]. However, it is presently unclear how 
these epigenetic changes affect the gene expression patterns 
of aging human tissues. 

We have previously used array-based methylation 
profiling of skin samples obtained from young and old 
volunteers and identified a defined set of genes that 
were hypermethylated in aged skin [16]. In the course 
of this study, two important features were identified 
that established the epidermis as an important tissue 
for further analysis: i) epidermis samples can be obtained 
by non-invasive procedures and are characterized by a 
very high degree of cell type homogeneity (>90% keratino- 
cytes), which greatly facilitates the identification of defined 
methylation patterns; and ii) the epidermis represents a 
tissue with high functional relevance for human health 
and disease and shows a well-known aging phenotype 
[17]. Furthermore, it has also been shown that the 
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DNMT1 DNA methyltransferase is functionally required 
for proper tissue homeostasis in the human epidermis [18]. 

Aging is a multifactorial process that results in a pro- 
gressive loss of regenerative capacity and tissue func- 
tionality. It has previously been suggested that these 
pathological changes are underpinned by a genome-wide 
loss of methylation marks [15]. We have now used 
transcriptome sequencing to identify genes that are dif- 
ferentially expressed in the young (18-24 years) and 
old (70-75 years) human epidermis. This revealed a 
highly specific set of 75 genes with strong functional 
relevance for human skin homeostasis. We also used 
whole-genome bisulfite sequencing to generate methy- 
lation maps from the same tissue samples. We find that 
the aging epidermis methylome is characterized by 
considerable stability and does not show any global 
methylation loss. However, local age-related methylation 
changes could be observed, and were enriched at weak 
or poised promoters and at enhancer regions, thus 
suggesting a functional relevance of epigenetic mechanisms 
for skin aging. 

Results 

The transcriptome of the aging human epidermis 

Primary human epidermis samples were obtained as 
suction blisters from the inner forearm of five female 
donors from two distinct age groups (18-24 years and 
70-75 years, see Additional file 1: Table SI for details), 
respectively. These two age groups provide a good 
representation of the skin aging phenotype in healthy 
adults [17]. After RNA purification, samples were pooled 
in equimolar ratios to allow library preparation for 
transcriptome sequencing. Pooling of samples was required 
to obtain sufficient amounts of mRNA for library prep- 
aration. Paired-end sequencing on an Illumina HiSeq 
2000 platform with read-lengths of 105 bases generated 
90 Gb of DNA sequence. After trimming to a maximal 
read length of 80 bases and a minimum base quality of 
a 30 Phred score, sequence reads were mapped to the 
GRCh37/hgl9 reference sequence using TopHat [19]. 
The resulting average strand-specific genome coverages 
were 500x for each sample (Table 1). 

Overall, the epidermis transcriptome was characterized 
by high expression levels of epidermis -specific genes, 

Table 1 Sequencing data 



such as various epidermis keratins and the epidermis 
structural genes Loricin and Filaggrin (Table 2). Genes 
that are predominantly expressed in cells of mesenchymal 
origin, such as Vimentin and Desmin, could not be 
detected or were only expressed at low levels (Table 2), 
which is in agreement with the substantial cell-type 
homogeneity of the samples used for this study. The 
expression patterns of these genes appeared highly 
similar in the young and old epidermis (Table 2), thus 
reflecting the stability of the core epigenetic program 
of differentiated human tissues. In agreement with this 
notion, gene expression levels for DNA methylation 
enzymes, such as DNMTs and TETs and their known 
cofactors DNMT3L and UHRF1, were low and also 
appeared very similar in the old and young samples 
(Table 2). 

To systematically identify genes that were differentially 
expressed in the young and the old sample, we used 
DESeq [20]. This revealed 75 genes with a statistically 
significant difference (Q <0.05, Additional file 2: Table S2), 
suggesting that the young and old epidermis transcriptomes 
were similar overall. This similarity was also confirmed 
by the relatively small effect size of age-related differential 
gene expression (Figure 1A). We also used an independent 
algorithm for differential expression analysis, Cuffdiff 2 
[21]. This identified an overall similar number of 107 
differentially expressed genes. Notably, out of the 75 
genes identified by DESeq, 68 were also identified by 
Cuffdiff 2 (Figure IB). This substantial overlap suggests 
that aging of the human epidermis affects the expression 
of a specific set of genes. 

Further analysis of differentially expressed genes iden- 
tified several factors that have previously been implied in 
skin homeostasis. For example, connective tissue growth 
factor (CTGF), which has been shown to be an import- 
ant factor in skin aging [22], was expressed at distinctly 
lower levels in the old epidermis (Figure 1C). We also 
used pathway analysis of differentially expressed genes 
to further characterize the functional effects of age- 
related differential gene expression. This revealed that 
genes involved in cell migration (P = 1.17 x 10" 16 ), cancer 
(P = 9.55 x 10" 13 ), dermatological diseases (P= 2.07 x 10~ n ) 
and cell proliferation (P = 1.4 x 10" 10 ), constituted the cat- 
egories with the most significant enrichment (Figure ID). 



Sample 


Age range 


Sequencing 


Number of read pairs 


Mapping efficiency 


Coverage 


Conversion rate 


Young (n = 5) 


18-24 


Transcriptome 


225,332,853 


89% 


500x 


n. d. 






Methylome 


344,185,333 


69% 


11. 3x 


99.84% 


Old (n = 5) 


70-75 


Transcriptome 


221,387,474 


88% 


500x 


n. d. 






Methylome 


334,873,828 


74% 


11. 9x 


99.88% 



Coverage indicates the average strand-specific genome or transcriptome coverage. Bisulfite conversion rates were determined by analyzing the conversion rates 
of co-purified mitochondrial genomes that were considered as unmethylated. Mitochondrial genome coverages were 1319x (young) and 1643x (old), respectively, 
n. d.: not determined. 
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Table 2 Expression levels of epidermis genes and DNA methylation factors 


Gene group 


Gene 


Expression (young) 


Expression (old) 


q value 


Epidermis keratins 


KRT1 


14719.2 


14954.2 


1.0 




KRT2 


7291.7 


6154.7 


1.0 




KRT5 


5860.6 


5247.7 


1.0 




KRT10 


26232.1 


23374.7 


1.0 




KRT14 


5665.0 


4430.4 


1.0 




KRT15 


382.1 


400.1 


1.0 


Epidermis structural factors 


LOR 


613.4 


427.3 


1.0 




FLG 


558.1 


501.3 


1.0 




FLG2 


322.5 


336.5 


1.0 


Fibroblast markers (dermis) 


VIM 


41.0 


54.0 


1.0 




DES 


<0.1 


<0.1 


1.0 


DNA methylation factors 


DNMT1 


9.5 


7.3 


1.0 




DNMT3A 


1.5 


1.2 


1.0 




DNMT3B 


1.4 


1.5 


1.0 




DNMT3L 


<0.1 


<0.1 


1.0 




UHRF1 


1.3 


1.2 


1.0 




TET1 


0.1 


0.1 


1.0 




TET2 


8.4 


6.9 


1.0 




TET3 


18.8 


14 


1.0 



Expression levels represent FPKM values [21]; levels >1 indicate expression; q values were calculated by DESeq [20]. 



In addition, functional annotation of differentially ex- 
pressed genes also showed a highly significant enrichment 
of genes involved in the development of the epidermis 
(P = 9.63 x 10~ 4 ), which also included downregulation of 
SPRR1A and B, KRT16 and KRT17 in old samples. In 
agreement with this finding, expression of genes involved 
in differentiation of keratinocytes was also altered in 
old samples, indicating an overall downregulation of 
differentiation (P = 1.43 x 10" 7 ). Altogether, our results 
thus suggest that deregulation of a relatively small set 
of genes could contribute to the phenotypic changes 
associated with human skin aging. 

To further analyze the patterns of differential gene 
expression, we experimentally validated the differential 
expression levels of defined genes in 18 individual epi- 
dermis samples from the two age groups (9 young and 
9 old). For these experiments, we selected 15 genes 
with a functional annotation related to skin homeostasis 
that had shown different degrees of age-related down- 
regulation in our RNA-seq analysis. The results showed 
significant (P <0.05, Mann-Whitney U test) differential 
gene expression for 8 out of 9 genes that had shown 
a > 2-fold expression change in our RNA-seq analysis 
(Figure 2). For genes with a fold change of <2, significant 
differential expression was only observed in 2 out of 6 
genes (Figure 2). These findings provide important 
confirmation for our transcriptome sequencing results 



and suggest that age-related gene deregulation occurs 
with a substantial degree of population homogeneity. 

The methylome of the human epidermis 

Having shown that aging is associated with the deregulation 
of a highly defined set of genes, we used whole-genome 
bisulfite sequencing to establish DNA methylation maps 
at single-base resolution. DNA was purified from the 
same epidermis samples that were used for transcriptome 
sequencing (Additional file 1: Table SI). Pooling of 
samples was necessary to achieve sufficient amounts 
of DNA for library preparation and has been previously 
used to reduce the effects of stochastic epigenetic variation 
[23]. Paired-end sequencing on an Illumina HiSeq 2000 
platform with read-lengths of 105 bases generated 137 Gb 
of DNA sequence. After trimming to a maximal read length 
of 80 bases and a minimum base quality of a 30 Phred 
score, sequence reads were mapped to the GRCh37/hgl9 
reference sequence using a mapping tool based on BSMAP 
2.0. The resulting average strand-specific genome coverage 
was 11.3x (young) and 11. 9x (old). We also determined 
the bisulfite conversion rate by analyzing mitochondrial 
sequences that were co-purified during the sample 
preparation and that we considered as unmethylated. These 
sequences showed a bisulfite conversion rate of >99.8% 
(Table 1), suggesting highly effective bisulfite treatment. 
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Figure 1 The transcriptome of the aging human epidermis. (A) MA plot showing fold change between young and old in relationship to 
the mean expression level. Red marks indicate genes with significant differential expression, as determined by DESeq analysis. (B) Venn diagram 
showing the overlap in differentially expressed gene sets, as determined by DESeq (green) and by Cuffdiff 2 (orange). (C) Connective Tissue 
Growth Factor (CTGF) is shown as a representative example for a differentially expressed gene. RNA-seq coverage is indicated for young (blue) 
and old (red) samples. (D) Ingenuity pathway analysis of differentially expressed genes. The plot shows the four most significantly enriched 
functional categories. 



Initial data analysis revealed that the human epidermis 
shares many basal features with published epigenomes 
from differentiated cultured human cell lines [3,4,24,25]. 
For example, the vast majority (>99.9%) of non-converted 
cytosines were found in a CpG dinucleotide context 
(Figure 3A), which is consistent with the overall de- 
amination efficiency and in agreement with the notion 
that non-CpG methylation is largely restricted to em- 
bryonic stem cells [24], Furthermore, methylation ratios 
of individual CpG dinucleotides revealed a characteristic 
bimodal distribution (Figure 3B). A major fraction of CpG 



dinucleotides (about 50%) showed complete methylation, 
as indicated by a methylation ratio of >0.95 (Figure 3B). 
Roughly 10% of the CpGs were completely unmethylated 
(methylation ratio <0.05), while 40% of the CpG dinu- 
cleotides showed partial methylation ratios between 
0.05 and 0.95 (Figure 3B). The average CpG methylation 
ratio was 0.74 (Figure 3C), which is again consistent with 
the CpG methylation ratios observed in other human 
datasets. Furthermore, average methylation ratios of 
promoter-associated CpGs were distinctly lower than 
the genome average, while gene bodies and intergenic 



Raddatz et al. Epigenetics & Chromatin 2013, 6:36 
http://www.epigeneticsandchromatin.eom/content/6/1/36 



Page 5 of 12 



gene 


Ioq 0 (fc) 


young 






old 




P-value 


CTGF 


3.0 












<0.01 


SPRR1A 


2.1 












0.01 


SPRR1B 


1.9 












0.01 


KRT17 


1.9 












<0.01 


KRT16 


1.5 












<0.01 


FZD10 


1.5 








■ 




0.02 


HA S3 


1.0 












<0.01 


PLIN2 


1.0 












0.09 


VCL 


1.0 


■ 


■ 


■ 






0.02 


RHPN2 


1.0 












0.11 


RGMA 


1.0 










0.57 


NPR3 


1.0 






■ H 


■ 




0.06 


VPS37B 


0.9 












0.01 


ERRFI1 


0.8 












<0.01 


PI4KB 


0.6 












0.66 




ACt 
















-2.5 








2.5 


Figure 2 Validation of differential age-related gene expression in individual tissue samples. qRT-PCR was performed on RNA from epidermal 


suction blisters of 18 healthy female volunteers (9 young and 9 old volunteers). The heatmap shows processed ACt values. Gene expression differences 


of individual samples are indicated relative to the gene-specific, age-independent average expression level over all samples (blue: lower ACt, red: higher 


ACt). Numbers in the left colum indicate fold-change differences in gene expression between young and old, as determined by transcriptome sequen- 


cing. P values were determined by a Mann-Whitney U test and indicate the significance for differential gene expression. The line between VCL and 


RHPN2 separates genes with a >2-fold change in gene expression from genes with a <2-fold change. 





regions showed higher methylation levels (Figure 3C), 
which is again similar to other published datasets. 

Overall, the methylation patterns of the young and old 
samples appeared very similar. This was evident not only 
by the average methylation ratios of individual genome 
compartments (Figure 3C), but also in comparisons of 
the global methylation landscapes (Figure 3D). A sliding 
window approach identified only 50 differentially meth- 
ylated windows of 100 kb (methylation difference >0.15), 
with an equal number of hypomethylated and hyper- 
methylated windows (Figure 3E). Similarly, a more local 
analysis with sliding windows of 5 CpGs did not reveal 
any directional changes in global methylation patterns 
(Figure 3F). Together, these findings strongly suggest 
that the global age-related methylation loss observed in 
T-cells [15] is not conserved in the epidermis. 

Identification and characterization of differentially 
methylated regions 

A visual inspection of the young and old methylation 
landscapes also indicated the presence of small clusters 
of differentially methylated CpG dinucleotides. To system- 
atically identify differentially methylated regions (DMRs), 
Fisher s exact test was used to determine the CpG dinucle- 
otides with a statistically significant (P <0.05) methylation 
difference. These differentially methylated CpGs (DMCs) 
were subsequently collapsed to identify regions of local, 
coordinated methylation changes. DMRs were defined as 
clusters of >8 DMCs with a distance of <50 bp between 
neighboring DMCs and a net region-wide methylation 



change of >8 DMCs. Only DMRs with an average sequen- 
cing coverage of >8 and methylation difference of >10% 
were used for further analysis. This identified 2,409 DMRs, 
of which 1,437 were more strongly methylated in the 
old sample, and 972 were more strongly methylated in 
the young sample. DMRs were comparably small (<150 bp) 
and associated with various gene regions. Notably, 1,156 
of these DMRs overlapped with the variably methylated 
regions that were recently identified through a compre- 
hensive analysis of 42 human methylomes [26]. This 
represents a robust (2.5-fold) enrichment over the genome 
average and suggests that age-related methylation changes 
affect epigenetic regulatory elements. 

To further characterize the DMRs, we used the avail- 
able ENCODE data [27,28] for normal human epidermal 
keratinocytes. The results showed that DMRs that were 
hypermethylated in the old sample were enriched 
for H3K27me3 and H3K4me3, the defining chromatin 
marks of poised promoters. In contrast, hypomethylated 
DMRs that were hypomethylated in the old sam- 
ple were enriched for H3K27Ac and also H3K4mel 
(Figure 4A), which represent established marks of en- 
hancers. We then used the available ChromHMM an- 
notation for normal human epidermal keratinocytes [29] 
to assign our DMR datasets to defined chromatin states 
(Additional file 3: Table S3). Subsequent data analysis 
showed a distinct enrichment of DMRs in promoters and 
enhancers (Figure 4B), which confirms and expands previ- 
ous observations on age-related methylation changes 
[13,14]. Furthermore, enhancer- associated DMRs showed 
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(See figure on previous page.) 

Figure 3 The methylome of the aging human epidermis. (A) Dinucleotide context of non-converted cytosine residues. (B) Methylation levels 
of individual CpG dinucleotides. Average methylation levels were determined for all covered CpG dinucleotides and then distributed into bins 
with increasing methylation ratios. Percentages indicate the fractions of unmethylated (light orange), partially methylated (orange) and completely 
methylated (dark orange) CpGs. (C) Average DNA methylation ratios of the genome (all), promoters, gene bodies and intergenic regions are shown for 
the young (blue) and old sample (red), respectively. (D) The methylation pattern of chromosome 8, shown in tracks of 100-kb windows. The 
blue line indicates the young sample, the red line indicates the old sample. (E) Density plot of average methylation ratios for 100-kb windows 
covering the entire genome. Numbers indicate the number of windows with a methylation difference >0.15 (dotted line). (F) Density plot of 
average DNA ratios for 5-CpG windows covering chromosome 4. 



a strong bias towards age-related hypomethylation 
(Figure 4B), which might reflect age-related activation 
of enhancers. Interestingly, active promoters showed 
no enrichment for DMRs (Figure 4B), which is again 
consistent with the comparably small effect size for age- 
related differential gene expression (Figure 1A). 

A visual inspection of DMRs provided further insight 
into their characteristic features. For example, the ERBB 
receptor feedback inhibitor 1 (ERRFI1) promoter region 
harbors a DMR in the shore region of the promoter- 
associated CpG island (Figure 4C). ERRFI1 is required 
for proper epidermal homeostasis [30] and was found 
to be hypermethylated and expressed at lower levels in 
old epidermis samples (Figure 4C). Another example 
for a DMR was identified in the low-density lipoprotein 
receptor (LDLR) gene region (Figure 4D), and was located 
in an annotated active enhancer element (Additional file 4: 
Figure SI). The DMR showed complex, but coordinated 
age-related methylation changes that were associated 
with lower expression levels of LDLR in the old sample 
(Figure 4D). Defects in the LDLR gene are the cause of 
familial hypercholesterolemia, which underlies the forma- 
tion of Xanthelasma [31], a dermatological lesion often 
found in the elderly population. These examples illustrate 
how age-related methylation changes can have relatively 
subtle, but significant expression changes on genes 
that are relevant for proper skin homeostasis. Further 
experiments will be required to determine the functional 
relevance of individual differentially methylated and 
expressed genes for the aging phenotype. 

Discussion 

Whole-genome bisulfite sequencing represents a powerful 
method to generate genome-wide methylation maps at 
single-base resolution [24,32]. This approach was re- 
cently used to characterize the methylomes of purified 
CD4 + T-cells from a newborn and a centenarian [15]. 
The results showed a pronounced destabilization of the 
aging epigenome, which was characterized by a widespread 
loss of methylation marks. Our analysis failed to detect 
any quantitative differences in the global methylation 
levels of young and old individuals. This can most likely 
be attributed to the different tissues used in both studies: 



while our analysis is based on primary epidermis samples, 
Heyn et al. [15] used purified CD4 + T-cells. Based on the 
specific features of T-cell differentiation and priming, it is 
conceivable that the epigenetic program of these cells is 
characterized by a particularly high degree of plasticity. 

On a global level, the transcriptomes and methylomes of 
the young and old epidermis appeared to be substantially 
similar. This is an important finding, because it illustrates 
the fundamental stability of the tissue-specific methylation 
landscape, which is required for the stable maintenance 
of cell type identity. Age-related methylation changes 
were limited to specific local alterations, which confirms 
and expands our previous observations [16]. With a size 
of 100-150 bp, age-related DMRs were considerably 
smaller than other known structures of the human methy- 
lome, such as DNA methylation valleys [3] and partially 
methylated domains [23-25,33,34], which extend over tens 
and hundreds of kilobases of DNA sequence, respectively. 
Further analysis will be required to understand the epigen- 
etic regulatory function(s) of these elements. 

Interestingly, our data also suggest that a significant 
fraction of the DMRs might represent enhancers that 
become aberrantly methylated during aging. The methyla- 
tion status of enhancer-associated CpG dinucleotides 
has recently been described to be closely associated with 
epigenetic gene deregulation in human cancers [35,36]. 
Our data identify similar elements at single-base resolution 
and suggest that the differential methylation of enhancers 
might be involved in age-related gene deregulation. 

Finally, our analysis identified several examples for 
hypermethylated DMRs in promoter regions. Promoter 
hypermethylation has been closely associated with gene 
silencing and plays an important role in the etiology of 
human tumors [37,38]. Furthermore, our data indicate 
an association of hypermethylated DMRs with bivalent 
chromatin structures. Similar results have been described 
in independent studies investigating age-related methyla- 
tion changes in other tissues [13,14], Bivalent chromatin 
modifications are a specific feature of stem cells and are 
not usually found in differentiated tissues like epidermis. 
The particular enrichment of age-related hypermethyla- 
tion in promoters that are annotated as "bivalent" might 
therefore reflect epigenetic changes in aging epidermal 
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Figure 4 (See legend on next page.) 
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Figure 4 Characterization of differentially methylated regions (DMRs). (A) Hypermethylated and hypomethylated DMRs show distinct 
chromatin states. DMRs that become hypermethylated in old epidermis (hyperDMRs) are enriched for H3K4me3 and H3K27me3, while DMRs that 
become hypomethylated in old epidermis (hypoDMRs) are enriched for H3K27ac. Numbers below the x-axis indicate the distance from the center 
of the DMR (dashed vertical line) in bp. (B) Enrichment of DMRs within defined genome segments. Bars indicate the ratio of the observed DMR 
frequency and the average frequency across the genome. Blue bars represent DMRs that are hypomethylated in the old epidermis, red bars represent 
DMRs that are hypermethylated in the old epidermis. (C) The ERBB receptor feedback inhibitor 1 (ERRFI1) promoter region harbors a representative 
DMR. ERRFI1 is required for proper epidermal homeostasis [30] and is expressed at lower levels in old epidermis samples (right panel, blue 
and red bars indicates expression in the young and old epidermis samples, respectively). Blue lines indicate methylation ratios in the young 
epidermis, red lines indicate methylation ratios in the old epidermis. The green bar indicates the position of the ERRFI1 promoter CpG island. (D) An 
annotated active enhancer element from the low-density lipoprotein receptor {LDLR) gene region harbors a DMR. Age-related methylation changes 
that were associated with lower expression levels of LDLR in the old sample (right panel, blue and red bars indicates expression in the young and old 
epidermis samples, respectively), which may promote the formation of Xanthelasma, a dermatological lesion often found in the elderly population. 



stem cells, and may underpin the decreased regenerative 
capacity of aging stem cells [39] . 

Conclusions 

Age-related DNA methylation changes have been described 
in various studies but their defining features and functional 
significance remain to be established. Our results show 
that the age-related genome-wide loss of DNA methy- 
lation observed in T-cells [15] is not conserved in the 
human epidermis, suggesting that the core developmental 
program of at least some human tissues is maintained 
through the aging process. Furthermore, our study 
represents the first to use transcriptome sequencing to 
explore the consequences of age-related epigenetic changes. 
Our results identified limited transcriptional changes 
that may underpin the aging phenotype. Finally, we 
identified several hundred highly localized elements 
with robust methylation differences between young and 
old. Interestingly, these elements were enriched for gene 
regulatory chromatin marks and many of them were 
associated with promoters and enhancers. Our study 
thus provides an important mechanistic framework for 
understanding age-related epigenetic deregulation. 

Methods 

Sample preparation for sequencing 

Epidermal suction blister samples were collected according 
to the current version of the Declaration of Helsinki 
and the guideline of the International Conference on 
Harmonization Good Clinical Practice (ICH GCP) was 
observed as applicable to a non-drug study. All volunteers 
provided written, informed consent. Suction blister 
samples were obtained at the study center of Beiersdorf 
AG and approved by the Beiersdorf AG Legal Review 
Board. Briefly, epidermis samples were detached from 
the forearms of 10 healthy female volunteers by applying 
a negative pressure of 150-250 mmHg for 2 h. Subse- 
quently, suction blister roofs were taken and immediately 
stored in liquid nitrogen. For nucleic acid isolation, 
suction blister samples were washed in DPBS (Cambrex, 
Verviers, Belgium) and homogenized using a TissueLyser 



(Retsch, Haan, Germany). DNA and RNA from suction 
blister samples were isolated using the QIAamp DNA 
Investigator Kit (Qiagen, Hilden, Germany) and RNeasy 
Fibrous Tissue Kit (Qiagen), respectively, according to 
the manufacturers instructions. The Poly(A)Purist MAG 
Kit (Ambion, Darmstadt, Germany) was used for mRNA 
selection. 

Sequencing 

Library preparation for bisulfite sequencing was performed 
as described previously [40]. Transcriptome sequencing 
libraries were prepared using the TruSeq RNA Sample 
Preparation Kit (Illumina, San Diego, USA), according to 
the manufacturers instructions. Paired-end sequencing 
was performed on an Illumina HiSeq system with read 
lengths of 105 base pairs and an average insert size of 
200 bp. 

Transcriptome mapping and quantification of 
differential expression 

RNA-seq reads were trimmed to a maximal length of 
80 bp and stretches of bases having a quality score <30 at 
the ends of the reads were removed. Reads were mapped 
using Tophat 2.0.6 [19]. As reference sequence for the 
transcriptome mapping we used the current assembly of 
the human genome (hgl9). Differential expression was 
quantified using DESeq 1.10.1 [20] applying the built-in 
procedures for library normalization and estimation of 
variance and with Cuffdiff 2.0 [21]. The resulting P values 
were subjected to multiple testing correction using built- 
in functions available in DESeq and Cuffdiff, respectively. 
Genes with a q value smaller than 0.05 were considered as 
differentially expressed. 

Pathway analysis 

Pathway analysis of differentially expressed genes was 
performed using IPA (Ingenuity Systems), using a P 
value cutoff of <0.05 on differential expression and a 
log fold-change of at least 0.263, corresponding to a 
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minimum expression change of 30% between young 
and old samples. 

qRT-PCR 

Total RNA was isolated from suction blister samples 
(n = 18) using the RNeasy Fibrous Tissue Kit (Qiagen) 
according to the manufacturer s instructions. After reverse 
transcription with the High Capacity cDNA Reverse Tran- 
scription Kit (Applied Biosystems, Darmstadt, Germany), 
samples were analyzed by TaqMan-PCR using the 
7900HT Fast-Real-Time PCR System (Applied Biosys- 
tems). The following assays were used, according to the 
manufacturers recommendations: CTGF (Hs01026927_gl), 
SPRR1A (Hs00954595_sl), SPRR1B (Hs00234164_ml), 
KRT16 (Hs00373910_gl), KRT17 (Hs01588578_ml), VCL 
(Hs00247826_ml), HAS3 (Hs00193436_ml), FDZ10 
(Hs00273077_sl), NPR3 (Hs00168558_ml), VPS37B 
(Hs00226582_ml), ERRFI1 (Hs00219060_ml), PLIN2 
(Hs00605340_ml), RHPN2 (Hs00369111_ml), RGMA 
(Hs00297192_ml), PI4KB (Hs01090927_ml). Data were 
analyzed utilizing the Sequence detector version 2.3 
software supplied with the 7900 Sequence Detector 
and RQ Manager 1.2. Quantification was achieved 
using the ACt method which indicates expression of 
the target gene relative to an endogenous reference 
(GAPDH; Hs99999905_ml). ACt values were averaged 
for all replicates of a gene/subject combination and for 
every gene the mean value over all subjects was sub- 
tracted, thus adjusting the average gene-specific ACt 
value to zero. The processed ACt values were visualized 
as heat maps. ACt values exceeding a threshold of 2.5 
were set to 2.5. 

Bisulfite mapping and methylation calling 

Reads were trimmed to a maximal length of 80 bp and 
stretches of bases having a quality score <30 at the ends 
of the reads were removed. Reads were mapped using 
BSMAP 2.02 [41]. As a reference sequence for the bisul- 
fite mapping we used the current assembly of the human 
genome (hgl9). Only reads mapping with both partners 
of the read pairs at the correct distance were used. The 
CpG-specificity was calculated by determining the number 
of cytosines called in all mapped reads at all non-CpG 
positions and dividing it by the number of all bases in 
all mapped reads at all non-CpG positions. Methylation 
ratios were determined using a Python script (methratio. 
py) distributed together with the BSMAP package. For 
both the forward and reverse strands, all cytosine bases 
in CG context were called independently. 

Identification and characterization of DMRs 

Fisher s exact test was used to identify 3,004,806 CpG 
dinucleotides with a statistically significant (P <0.05) 
difference in their methylation ratios between the young 



and the old datasets. Differentially methylated CpGs 
were subsequently collapsed to identify regions of local, 
coordinated methylation change. DMRs were defined 
as clusters of >8 DMCs with a distance of <50 bp between 
neighboring DMCs and a net region-wide methylation 
change of >8 DMCs. Only DMRs with an average 
sequencing coverage of >8 and methylation difference 
of >10% were used for further analysis. ChlP-seq peaks 
of histone modification marks in NHEK cells were 
downloaded from the ENCODE website [29] and averaged 
in 50 bp windows surrounding the DMR sites. Loess 
smoothing with a span of 10% was applied. For the en- 
richment analysis of DMRs in defined chromatin states, 
the observed frequency (%) of DMRs centered in given 
states (ChromHMM annotation data in NHEK cells, 
ENCODE website) was divided by the genomic frequency 
of the state. 

Data access 

Sequencing data have been deposited in the GEO database 
under the accession number GSE46487. 

Additional files 

f N, 

Additional file 1: Table SI. Samples used for sequencing. The table 
shows a complete overview of all epidermis samples used for methylome 
and transcriptome sequencing. 

Additional file 2: Table S2. List of differentially expressed genes. 
Differentially expressed genes were identified by DESeq using built-in 
procedures for library normalization and estimation of variance. baseMean is 
the average of baseMean (young) and baseMean (old) and indicates 
the number of fragments per gene after library normalization, P values 
were subjected to multiple testing correction using built-in functions 
available in DESeq, leading to adjusted P values ("padj"). Genes with an 
adjusted P value <0.05 were considered as differentially expressed and 
are included in this table. 

Additional file 3: Table S3. Association of DMRs with ChromHMM 
segments. The table shows the results obtained from ChromHMM 
segmentation of the human genome sequence, using ENCODE data for 
normal human keratinocytes. 

Additional file 4: Figure SI. Epigenomic analysis of the LDLR gene 
region. UCSC Genome Browser tracks for DNA methylation and various 
histone marks, based on ENCODE data for normal human keratinocytes. 
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